package edu.gmu.atelier;

import java.util.ArrayList;
import java.util.List;

/**
 * Performs the Simplex Method to determine optimal solution to a LP,
 * linear program.
 * 
 * @author James H. Pope
 */
public class SimplexSolver implements Solver
{
    private Matrix A = null;
    private Vector b = null;
    private int[] B  = null;
    
    public SimplexSolver( Matrix A, int[] B, Vector b )
    {
        this.A = A;
        this.B = B;
        this.b = b;
    }
    
    public Matrix getA()
    {
        return this.A;
    }
    
    
    private static Vector solve( Matrix A, int[] B, Vector c, Vector b )
    {
        int iteration = 0;
        while(iteration < 100)
        {
            VectorFactory f = A.getFactory();
            Matrix Ab = A.extract(B);
            Vector bb = b.extract(B);
            System.out.println(""+iteration+"===============================" );
            System.out.println("Ab: "+ Ab );
            System.out.println("bb: "+ bb );

            /*
             * Ab x* = bb solve for x*
             * 
             * x* is the current basic feasible solution
             */
            LUFactor abf = new LUFactor(Ab);
            Vector xStar = abf.solve(bb);
            System.out.println("x*: "+ xStar );

            /*
             * Objective function value produced at x*
             */
            double objValue = c.dot(xStar).doubleValue();
            System.out.println("cTx: "+ objValue );

            /*
             *  lbT Ab=cT
             * Ab^T ib=c
             */
            LUFactor abTf = new LUFactor(Ab.transpose());
            Vector lb = abTf.solve(c);
            System.out.println("lb: "+ lb );

            /*
             * Which index leaves the basis?  The minimum one.
             */
            double minValue = Double.MAX_VALUE;
            int iExit = -1;
            boolean isNegative = false;
            for (int i = 0; i < lb.getLength(); i++)
            {
                double value = lb.getElement(i).doubleValue();
                if( value < minValue )
                {
                    minValue = value;
                    iExit = i;
                }
                if( value < 0 )
                {
                    isNegative = true;
                }
            }
            System.out.println("Leaving index i: "+ iExit );
            // Check for optimality
            if( isNegative == false )
            {
                return xStar;
            }

            /*
             * Now determine which direction d to head from x*
             * This is Ab d = { 0,-1, 0}
             *            d = { 0,-1, 0} Ab.inv()
             * This is inefficient compared to ABi.iExit * -1
             */
            Vector abd = f.makeVector( bb.getLength() );
            abd.setElement(iExit, f.makeOne().negate());
            Vector d = abf.solve(abd);
            System.out.println("Moving in direction d: "+ d );

            /*
             * Now have to do some work to determine index to enter basis
             */
            /*
             * We are determining the value alpha * d of how 
             * far we should move in direction d before hitting 
             * another adjacent constraint.
             * This is when d * normal inequality is positive
             * (going similar direction), computed by Ad
             */
            Vector ad = A.mult(d);
            System.out.println("ad: "+ ad );

            /*
             * Determine index set K, i.e. those that are stritcly positive
             * Note that if K is empty, then the LP is unbounded
             */

            List<Integer> K = new ArrayList<Integer>();
            for (int i = 0; i < ad.getLength(); i++)
            {
                if( ad.getElement(i).doubleValue() > 0 )
                {
                    K.add(i);
                }
            }
            System.out.println(" K: "+ K );

            if( K.isEmpty() )
            {
                System.out.println("LP is unbounded");
                return null;
                // might also throw exception
            }



            /*
             * Which one enters from the index set K?
             * 
             * min{ bk - akT x* / akT d }
             * 
             * Said another way, the one with the smallest positive
             * (i.e. going in obj function direction) value, which 
             * represents the first inequality to current vertex.
             */
            // Crappy inefficient way
            Vector bAx = b.sub( A.mult(xStar) );
            double minKValue = Double.MAX_VALUE;
            int jEnter = -1;
            for (int i = 0; i < K.size(); i++)
            {
                int k = K.get(i);
                double kValue = bAx.getElement(k).doubleValue() /
                                ad.getElement(k).doubleValue();
                System.out.println("kValue = "+kValue);
                if( kValue < minKValue )
                {
                    minValue = kValue;
                    jEnter= k;
                }
            }
            System.out.println("Entering index j: "+jEnter);

            // Swap out i and j
            B[iExit] = jEnter;
            System.out.println("Updated B: "+ print(B) );
            
            /*
             * Note: Python version updates Ab and Ab^-1 in an efficient
             * way O(n^2) instead of recomputing.  Ab no problem, simply
             * a row replacecment.  Ab^-1 much harder and we currently
             * do not support matrix inverses (always rearrange to use
             * LU factorization and then solve).
             */
            
            iteration++;
        }
        // There is a proof that the algorithm ends
        return null;
    }
    
    public static final int[] makeRange( int start, int length )
    {
        int[] B = new int[length];
        for (int i = 0; i < B.length; i++)
        {
            B[i] = start+i;
        }
        return B;
    }
    
    public static final String print( int[] a )
    {
        StringBuilder b = new StringBuilder("[");
        for (int i = 0; i < a.length; i++)
        {
            b.append(a[i]+",");
        }
        b.append("]");
        return b.toString();
    }
    
    //------------------------------------------------------------------------//
    // Interface methods
    //------------------------------------------------------------------------//
    public Vector solve( Vector c )
    {
        /*
         * Preconditions:
         * 1) A is linearly independent (rank(A)=n)
         * 2) Non-degenerate, which means basic feasible solutions and
         * basis "subsets" are one-to-one (degenerate can have many 
         * basis subsets for a single vertex, i.e. basic feasible
         * solution).
         */
        //int[] B = new int[] {0,1,2};
        //int[] B = new int[] {0,5,2};
        Vector xOpt = solve(A,B,c,b);
        return xOpt;
    }
    
}